// @HEADER
//*********************************************************************//
//  SiYuan: A numerical PDE solver                                     //
//  Copyright (2022) YUAN Xi                                           //
//  This Software is released under the BSD 2-Clause license detailed  //
//  in the file "LICENSE" in the top-level SiYuan directory            //
//*********************************************************************//
// @HEADER

#include "InitialCondition.hpp"

#include "Thyra_VectorStdOps.hpp"

#include "Panzer_GlobalEvaluationDataContainer.hpp"
#include "Panzer_STK_SetupUtilities.hpp"
#include "Panzer_ThyraObjContainer.hpp"

template<typename ScalarT>
void SiYuan::InitialConditions(Teuchos::ParameterList& pl, 
    Thyra::ModelEvaluatorDefaultBase<ScalarT> & model,
    Teuchos::RCP<panzer_stk::STK_Interface>& mesh,
	Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> >& lof,
    Teuchos::RCP<panzer::PhysicsBlock>& physicsBlock)
{
	PHX::FieldManager<panzer::Traits> fm;
	Teuchos::RCP<panzer::LinearObjContainer> loc = lof->buildGhostedLinearObjContainer();

    /*Teuchos::RCP<PHX::Evaluator<panzer::Traits> > evaluator = lof->buildGather<panzer::Traits::Residual>(pl);
	fm.registerEvaluator<panzer::Traits::Residual>(evaluator);
	fm.requireField<panzer::Traits::Residual>(*evaluator->evaluatedFields()[0]);*/

    Thyra::ModelEvaluatorBase::InArgs<double> nomValues = model.getNominalValues();
    Teuchos::RCP<Thyra::VectorBase<double> > x_vec = Teuchos::rcp_const_cast<Thyra::VectorBase<double> >(nomValues.get_x());
    Thyra::seed_randomize<double>(123456u);
    Thyra::randomize(0.0,1.0,x_vec.ptr());
    Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(loc)->set_x_th(x_vec);

    std::vector<bool> active_evaluation_types(Sacado::mpl::size<panzer::Traits::EvalTypes>::value,false);
    int residual_index = Sacado::mpl::find<panzer::Traits::EvalTypes,panzer::Traits::Residual>::value;
    active_evaluation_types[residual_index] = true;
    physicsBlock->setActiveEvaluationTypes(active_evaluation_types);
    Teuchos::ParameterList user_data;
    physicsBlock->buildAndRegisterGatherAndOrientationEvaluators(fm, *lof, user_data);
	
	panzer::Traits::SD sd;
    fm.postRegistrationSetup(sd);

    panzer::Traits::PED ped;
    ped.gedc->addDataObject("Initial Condition Gather Container",loc);
//    ped.gedc->addDataObject("Residual Scatter Container",loc);
    fm.preEvaluate<panzer::Traits::Residual>(ped);

    // run tests
    /////////////////////////////////////////////////////////////
    Teuchos::RCP<std::vector<panzer::Workset> > work_sets = 
        panzer_stk::buildWorksets(*mesh,physicsBlock->elementBlockID(), physicsBlock->getWorksetNeeds()); 

    panzer::Workset & workset = (*work_sets)[0];
    workset.alpha = 0.0;
    workset.beta = 1.0; // derivatives multiplied by 2
    workset.time = 0.0;
    workset.evaluate_transient_terms = false;

    fm.evaluateFields<panzer::Traits::Residual>(workset);
}
